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Brown [Phys. Rev. D 79, 104029 (2009)] has recently introduced a covariant formulation of the 
BSSN equations which is well suited for curvilinear coordinate systems. This is particularly desirable 
as many astrophysical phenomena are symmetric with respect to the rotation axis or are such that 
curvilinear coordinates adapt better to their geometry. However, the singularities associated with 
such coordinate systems are known to lead to numerical instabilities unless special care is taken 
(e.g., regularization at the origin). Cordero-Carrion will present a rigorous derivation of partially 
implicit Runge-Kutta methods in forthcoming papers, with the aim of treating numerically the stiff 
source terms in wave-like equations that may appear as a result of the choice of the coordinate 
system. We have developed a numerical code solving the BSSN equations in spherical symmetry 
and the general relativistic hydrodynamic equations written in flux-conservative form. A key feature 
of the code is that it uses a second-order partially implicit Runge-Kutta method to integrate the 
evolution equations. We perform and discuss a number of tests to assess the accuracy and expected 
convergence of the code, namely a pure gauge wave, the evolution of a single black hole, the evolution 
of a spherical relativistic star in equilibrium, and the gravitational collapse of a spherical relativistic 
star leading to the formation of a black hole. We obtain stable evolutions of regular spacetimes 
without the need for any regularization algorithm at the origin. 
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I. INTRODUCTION 

The 3+1 formulation of Einstein equations originally 
proposed by Nakamura [l[ and subsequently modified by 
Shibata-Nakamura 2] and Baumgarte-Shapiro Q, which 
is usually known as the BSSN formulation, has become 
the most widespread used formulation in the numerical 
relativity community. This is due to its stability prop- 
erties, and to the developments associated with gauge 
conditions and the puncture method which have proved 
essential to perform accurate and long-term stable evo- 
lutions of spacetimes containing black holes (BHs) [HIH- 

The main drawback of the BSSN formulation in its 
original form resides in the fact that it is particularly 
tuned for Cartesian coordinates, since this involves dy- 
namical fields which are not true tensors and assumes 
that the determinant of the conformal metric is equal to 
one. Brown @ addressed this issue and introduced a co- 
variant formulation of the BSSN equations which is well 
suited for curvilinear coordinate systems. This is par- 
ticularly desirable as many astrophysical phenomena are 
symmetric with respect to the rotation axis (e.g., accre- 
tion disks) or are such that spherical coordinates adapt 
better to their geometry (e.g., gravitational collapse). 

However, the singularities associated with the curvilin- 
ear coordinate systems are a known source of numerical 
problems. For instance, one problem arises because of the 
presence of terms in the evolution equations that behave 
like 1/r near the origin r = 0. Although on the analyt- 
ical level the regularity of the metric ensures that these 
terms cancel exactly, on the numerical level this is not 
necessarily the case, and special care should be taken in 
order to avoid numerical instabilities. A similar problem 



appears also near the axis of symmetry in axisymmetric 
systems if curvilinear coordinate systems are used. 

Several methods have been proposed to handle the is- 
sue of regularity in curvilinear coordinates. One possible 
approach is to rely on a specific gauge choice (i.e., the 
polar/ areal gauge) @, [sfj but it has the obvious limita- 
tion of restricting the gauge freedom which is one of the 
main ingredients for successful evolutions with the BSSN 
formulation. An alternative method is to apply a regu- 
larization procedure. One such regularization methods, 
presented by [9[ , enforces both the appropriate parity reg- 
ularity conditions and local flatness in order to achieve 
the desired regularity of the evolution equations. Such 
method has the advantage that it allows a more generic 
gauge choice, and has been explored by [ToT - fT^ I who have 
performed several numerical simulations of regular space- 
times in spherical and axial symmetry. In particular, 
in [l2j , the authors applied a regularization algorithm to 
the BSSN equations in spherical symmetry. A disadvan- 
tage of such a regularization algorithm is that it is not 
easy to implement numerically both conditions simulta- 
neously, and it requires the introduction of auxiliary vari- 
ables as well as finding their evolution equations. This 
is an obstacle if one wants to perform 3D simulations of 
regular spacetimes with spherical coordinates. 

Therefore, one would ideally like to use a numerical 
scheme that is able to integrate in time a system of equa- 
tions like the BSSN, in curvilinear coordinates (with or 
without symmetries), without the burden of regulariza- 
tion in order to achieve the desired stability and robust- 
ness. Implicit or partially implicit methods are used to 
deal with systems of equations that require a special nu- 
merical treatment in order to achieve stable evolutions. 
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The origin of the numerical instabilities may be diverse. 
Stiff source terms in the equations can lead to the devel- 
opment of numerical instabilities, and with some choices 
of the coordinate system, source terms may introduce fac- 
tors which can be numerically interpreted as stiff terms 
(e.g., 1/r factors due to spherical coordinates close to 
r = even when regular data is evolved). Recently, par- 
tially implicit Runge-Kutta (PIRK) methods for wave- 
like equations in spherical coordinates have been suc- 
cessfully applied [l3| to the hyperbolic part of Einstein 
equations in the Fully Constrained Formulation [bH ]. 

The first steps through the rig orous derivation of the 
PIRK methods will appear in [15[ and a detailed descrip- 
tion of the methods and their properties will be derived in 
a forthcoming paper [16j . Motivated by these results, we 
have developed a numerical code solving the BSSN equa- 
tions in spherical symmetry and the general relativis- 
tic hydrodynamics equations written in flux-conservative 
form [rjj. The code uses a second-order PIRK method 
to integrate the evolution equations in time, and we do 
not apply any regularization scheme at the origin. This 
approach has the additional advantages that it imposes 
no restriction at all on the gauge choice (one can there- 
fore use the moving puncture gauge) and no special care 
should be taking in the transition between a regular 
spacetime and that containing a singularity as it hap- 
pens in the gravitational collapse of a star to a BH. 

The paper is organized as follows. The formulation 
of Einstein equations, including the implementation of 
the puncture approach and gauge conditions, along with 
the formulation of the general relativistic hydrodynamic 
equations is briefly presented in Sec. II. Sec. Ill gives a 
short description of the PIRK method used, while Sec. IV 
describes the numerical implementation. Sec. V discusses 
numerical simulations of a pure gauge wave, the evolu- 
tion of a single BH, the evolution of spherical relativistic 
stars in equilibrium, and the gravitational collapse of a 
spherical relativistic star leading to the formation of a 
BH. A summary of our conclusions is given in Sec. VII. 
We use units in which c = G = M Q = 1. Greek indices 
run from to 3, Latin indices from 1 to 3, and we adopt 
the standard convention for the summation over repeated 
indices. 



II. BASIC EQUATIONS 

We next give a brief overview of the formulation for 
the system of Einstein and hydrodynamic equations as it 
has been implemented in the code. 



A. BSSN equations in spherical symmetry 

A reformulation of the ADM system, the BSSN for- 
mulation |1H3J, has been implemented to solve Einstein 
equations. In particular, we solve the BSSN equations in 



the special case of spherical symmetry. We refer to [12J 
for a detailed description of the equations. 

Under this symmetry condition the spatial line clement 
is written as 

dl 2 = e 4x [a(r,t)dr 2 +r 2 b{r,t)dn 2 ], (2.1) 

where dil 2 is the solid angle element, d£l 2 = d9 2 + 
sin 2 9dtp 2 , a(r, t) and b(r, t) are the metric functions, and 
X is the conformal factor defined as 

X = ^M7/7), (2-2) 

where 7 is the determinant of the conformal metric. The 
conformal metric relates to the physical one by 

% = e" 4 >%-- (2.3) 

Initially, the determinant of the conformal metric ful- 
fills the condition that it equals the determinant of the 
flat metric in spherical coordinates jij (i.e. j(t = 0) = 
7 = r 4 sin 2 0). Moreover, we follow the so called "La- 
grangian" condition dt'J = (i.e. choosing a = 1 in 
Eqs. (2.4), (2.6), (2.7), (2.16) and (2.17)). The evolution 
equation for the conformal factor takes the form 

d tX = p r d rX + cjV m l3 m -l;aK, (2.4) 

6 

K being the trace of the extrinsic curvature, a the lapse 
function, and 

the divergence of the shift vector f3 l . The evolution equa- 
tions for the conformal metric components are: 

d t a = (3 r d r a + 2ad r p r - ^aaV m (3 m - 2aaA a , (2.6) 

d t b = f3 r d r b + 2b— - -abV m f3 m - 2abA b , (2.7) 
r 3 

where Aij is the traceless part of the conformal extrinsic 
curvature, and 

A a = A r r , A b = Ag. (2.8) 

Note that as A^ is traceless A a + 2A b — 0. The evo- 
lution equation for K is: 

d t K = f3 r d r K - V 2 a + a(A 2 a + 2A 2 + \k 2 ) 

+ 4Tra(E + S a + 2S b ), (2.9) 

with the matter source terms measured by the Eulerian 
observers given by 

E = n^n v T»\ 
ji = -li^n v T^ v , 
Sij = H^juT^, (2.10) 
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T^ v being the stress-energy tensor for a perfect fluid, 
which is written as a function of the rest-mass density 
p, the specific enthalpy h, the pressure P and the fluid 
4- velocity u M , 

= phu"u v + Pg^, (2.11) 

and 

S a = , Sb = Sg. (2-12) 

The Laplacian of the lapse function with respect to the 
physical metric is given by 



In addition to the evolution equations there are con- 
straint equations, the Hamiltonian and the momentum 
constraints, which are only used as diagnostics of the ac- 
curacy of the numerical evolutions: 

n = R- {A 2 a + 2A 2 b ) + \K 2 -WttE = 0, (2.18) 
M r = d r A a - 2 -dvK + 6A a d rX 

+ (A a - A b ) fl + M\- SttjV = 0. (2.19) 



V 2 a 



, I d r a d r b 2 

d r a -d r a — - 2d r x 

2a b r 



(2.13) 



Next, the evolution equation for the independent com- 
ponent of the traceless part of the conformal extrinsic 
curvature, A a , is given by 



V r V r a - l^a) +a(R r r - ^R 



d t A a = (3 r d r A a 

+ aKA a — l6ira(S a — Sb) 



(2.14) 



where i?£ is the mixed radial component of the Ricci 
tensor, R its trace, and V r V r o; is written as 



V r V r a 



1 



o4 X 



>>;n ■ 20, \ 



(2.15) 



Finally, the evolution equation for A r , the radial com- 
ponent of the additional BSSN variables A* 
with A a bc = f a bc - f g c , is given by 



d t A r 



r d r A r - k r d T p r + -d 2 f3 r + %d r 
a b 



3 V a 



-(A a d r a + ad r A a ) 
a 



+ 2a A a A 



rb 



[A a - A b 



+ 



drA a - TprK + QA a 8 r X 



, , ,'2 d r b . 

+ (A a -A h )[- + -^-\- %TT 3r 



(2.16) 



where we take £ = 2. 

Note that in the simulations shown in Sec. V we have 
evolved the quantity X = e~ 2x instead of the conformal 
factor x (although similar conclusions can be drawn if the 
conformal factor x is used instead). We replace Eq. (|2.4[) 
by the following evolution equation for X: 

d t X = p r d r X - ^X(aK - aV m f). (2.17) 



1. Gauge choices 

In addition to the BSSN spacetime variables, there are 
two more variables left undetermined, the lapse, a, and 
the shift vector, f3 l . The code can handle arbitrary gauge 
conditions, however unless otherwise indicated, we use 
the so called "non-advective 1+log" condition [18{ for the 
lapse, and a variation of the "Gamma-driver" condition 
for the shift vector (HEH. 

The form of this slicing condition is expressed as 



d t a = -2aK. 



(2.20) 



For the radial component of the shift vector, we choose 
the Gamma-driver condition, which is written as 



d t B r = -d t A r 
i 4 1 



(2.21) 

d t (3 r = B\ (2.22) 
where the auxiliary variable B r is introduced. 

B. Formulation of the hydrodynamic equations 

The general relativistic hydrodynamic equations, ex- 
pressed through the conservation equations for the stress- 
energy tensor T^ iV and the continuity equation are: 



V M T^ = , V p (pi*") = 0. 



(2.23) 



Following 17], the general relativistic hydrodynamic 
equations are written in a conservative form in spherical 
coordinates. The following definitions for the hydrody- 
namic variables are used: 



fir 

a 



W = au\ 



(2.24) 



(2.25) 



where W is the Lorentz factor. By defining the vector of 
unknowns, U, as 



U= y/j(D,Sr,T), 



(2.26) 
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where the conserved quantities are 

D = pW, 
S r = phW 2 v r , 
t = phW 2 -P-D, 

and fluxes, F r , as 

F r = 



= ^—g[D{v r -p r /a), 
S r (v r - (3 r /a) + P, 
r(v r - /3 r /a) + Pv r ] , 



(2.27) 
(2.28) 
(2.29) 



(2.30) 



the set of hydrodynamic equations (|2 .23[) can be written 
in conservative form as 



flfcU + d r F r = S, 
where S is the vector of sources given by 
1 



0,T 00 ( -(^fd rlrr - a d r a 



+T 0r /3 r d r ~/ rr + T?d r f3 r + -T rr d rlrri 
(T 00 /3 r + T 0r^pr Krr _ + T rr Rrr 



(2.31) 



(2.32) 



To close the system of equations, we choose the T-law 
equation of state given by 

P=(T-l)pe, (2.33) 

where e is the specific internal energy. 

III. PIRK METHODS 

Let us consider the following system of PDEs, 



u t = Ci(u,v), 

v t = £ 2 {u) +£ 3 (u,v), 



(3.1) 



Ci, C 2 and £3 being general non-linear differential op- 
erators. Let us denote by L\, L 2 and L3 their discrete 
operators, respectively. L\ and L3 will be treated in an 
explicit way, whereas the L 2 operator will be considered 
to contain the unstable terms and, therefore, treated par- 
tially implicitly. 

We use a Runge-Kutta (RK) method to update in 
time the previous system (|3.ip . Each stage of the PIRK 
method consists of two steps: i) the variable u is evolved 
explicitly; ii) the variable v is evolved taking into account 
the updated value of u for the evaluation of the L 2 oper- 
ator. This strategy implies that the computational costs 
of the methods are comparable to those of the explicit 
ones. The resulting numerical schemes do not need any 
analytical or numerical inversion, but they are able to 
provide stable evolutions due to their partially implicit 
component. 



For the numerical simulations shown in the paper, we 
use the second-order PIRK scheme, which follows as: 



AtL 1 {u , \v n ), 



At 



±L 2 (u n ) + ±L 2 (u^)+L 3 (u n ,v n ) 



(3.2) 



,71+1 



,,7i+l 



i B + t»W+Atii(u( 1 ),t; (1) ) 



At 



v n + — [L 2 (u n ) + L 2 (u n+L ) 

" +L 3 (u n ,v n ) + L 3 (uW,vW) 

(3.3) 

In the first stage, u is evolved explicitly; the updated 
value vP' is used in the evaluation of the L 2 operator for 
the computation of Once all the values of the first 
stage are obtained, we proceed to the final one. Again, 
u is evolved explicitly (using the values of the variables 
of the previous time-step and previous stage), and the 
updated value u n+1 is used in the evaluation of the L 2 
operator for the computation of v n+1 . 

This scheme is applied to the hydrodynamic and BSSN 
evolution equations. We include all the problematic 
terms appearing in the sources of the equations in the 
L 2 operator. Firstly, the hydrodynamic conserved quan- 
tities, the conformal metric components, a and 6, the 
conformal factor, x, or the quantity X (function of the 
conformal factor), the lapse function, a, and the radial 
component of the shift, f3 r , are evolved explicitly (as u 
is evolved in the previous PIRK scheme); secondly, the 
traceless part of the extrinsic curvature, A a , and the trace 
of the extrinsic curvature, K , are evolved partially im- 
plicitly, using updated values of a, a and b; then, the 
quantity A r is evolved partially implicitly, using the up- 
dated values of a, a, b, /3 r , conformal factor, A a and 
K; finally, B r is evolved partially implicitly, using the 
updated values of A r . Matter source terms are always 
included in the explicitly treated parts. In Appendix [X[ 
we give the exact form of the source terms included in 
each operator. 

The PIRK methods will be further described in a forth- 
coming paper (l6j and derived up to third-order in At 
(time-step), in such a way that the number of stages is 
minimized. These methods are based on stability prop- 
erties for both the explicit and implicit parts, recover- 
ing the optimal SSP explicit RK methods [2(| when the 
L 2 operator is neglected, i.e., partially implicitly treated 
parts are not taken into account. 



IV. IMPLEMENTATION 

A. Numerics 

Derivatives in the spacetime evolution equations are 
calculated using a fourth-order centered finite difference 
approximation in a uniform grid except for the advection 
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terms (terms formally like j3 r d r u) 1 for which an upwind 
scheme is used. We also use fourth-order Kreiss-Oliger 
dissipation [2l[ to avoid high frequency noise appearing 
near the outer boundary. 

We use a second-order slope limiter reconstruction 
scheme (MC limiter) to obtain the left and right states of 
the primitive variables at each cell interface, and a HLLE 
approximate Riemann solver (22l |23| . 






\f\ 


1 

K, 1=5 
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B. Boundary Conditions 

The computational domain is defined as < r < L, 
where L refers to the location of the outer boundary. 
We used a cell-centered grid to avoid that the location 
of the puncture at the origin coincides with a grid point 
in simulations involving a BH. At the origin we impose 
the conditions derived from the assumption of spherical 
symmetry. At the outer boundary we impose radiative 
boundary conditions [l9j for the spacetime variables ex- 
pressed as 
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FIG. 1: Trace of the extrinsic curvature, K, for a pure gauge 
pulse as a function of the radius at four different times. 
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Off = -Vdrf - -(f - f ), 

r 



(4.1) 



where /o is the background solution of the field and v is 
the wave speed. 



C. Atmosphere treatment 

An important ingredient in numerical simulations 
based on finite difference schemes to solve the hydrody- 
namic equations is the treatment of vacuum regions. The 
standard approach is to add an atmosphere of very low 
density filling these regions [24| . We follow this approach 
and treat the atmosphere as a perfect fluid with a rest- 
mass density several orders of magnitude smaller than 
that of the bulk matter. The hydrodynamic equations 
are solved in the atmosphere region as in the region of 
the bulk matter. If the rest-mass density p or specific 
internal energy e fall below the value set for the atmo- 
sphere, these values are reset to have the atmosphere 
value of the primitive variables. 



V. NUMERICAL RESULTS 

A. Pure gauge dynamics 

We first consider the propagation of a pure gauge pulse 
using the same initial parameters as in |12j |. The main 
difference with respect to [l2j is that we do not regularize 
the origin and rely only on the PIRK scheme to achieve 
a stable numerical simulation. The initial data are given 



by 



X 


= 0, 




a 


= b = 


1, 


A a 


= A b 


= K = 


A r 


= 0, 




a 


= 1 + 


a r 2 




1 + r 2 



<r-r ) 2 + e -(r+r ) 2 



(5.1) 

(5.2) 

(5.3) 
(5.4) 

(5.5) 



with ao = 0.01 and ro = 5. We evolve these initial data 
with a grid resolution of Ar = 0.1 and At — 0.5Ar. We 
use zero shift and harmonic slicing, 



d±a = 



-a 2 K. 



(5.6) 



In Fig. [TJ we show the trace of the extrinsic curvature, 
K, as a function of the radius at four different times 
(t = 0, 5, 10, 15). The initial pulse separates in two pulses 
propagating in opposite directions. The snapshots of the 
evolution of the trace of the extrinsic curvature show that 
the evolution remains well behaved everywhere in the 
computational grid. We note that at t = 5 the value of 
K reaches a value of 0.1 at the origin, but later returns 
to zero when the pulse moves outwards as shown by [12j . 

In Fig. [3J we plot the Hamiltonian constraint 
(Eq.flnHJ}) at four different times (t = 0,5,10,15). Al- 
though the largest violation of the Hamiltonian con- 
straint occurs close to the origin (and is of the order 
of 10~ 3 for the time frames shown in Fig. [5]), we find 
that it remains well behaved and there is no sign of any 
numerical instability despite the fact the initial data are 
regular at the origin and we do not impose any regularity 
conditions there. 
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FIG. 2: Hamiltonian constraint for a pure gauge pulse as a 
function of the radius at four different times. 



In order to asses the convergence of the code, we have 
performed three simulations with resolutions Ar = 0.1, 
Ar = 0.05 and Ar = 0.025. The Hamiltonian constraint 
violations rescaled by the factors corresponding to second 
order convergence at t = 10 are plotted in Fig. [3l All 
three lines overlap indicating that the code achieves the 
second-order convergence expected for the PIRK scheme 
used. 



B. Schwarzschild black hole 

The Schwarzschild metric in isotropic coordinates is 
used as initial data to test the ability of the code to evolve 
BH spacetimes within the moving puncture approach. 
The initial data are such that the 3-metric is written as 



dl 2 = ijj 4 (dr 2 +r 2 dn 2 ), 



(5.7) 



where the conformal factor is tp = (1 + M/2r), M being 
the mass of the BH, which we set as M = 1. Here r is 
the isotropic radius. Initially the extrinsic curvature is 
Kn = 0. 

We evolve the single stationary puncture initial data 
with a precollapsed lapse and initially vanishing shift vec- 
tor. We use the gauge conditions given by Eqs. (|2.20l) - 
(|2~22|) . with a resolution Ar = 0.05, At = 0.5Ar and 
N r = 30000 grid points to place the outer boundary suf- 
ficiently far way from the puncture so that errors from 
the boundary do not affect the evolution. 

As pointed out by [25|, [26[, the numerical slices of a 
Schwarzschild BH spacetime with these gauge conditions 
reach a stationary state after t ~ 20. This is shown in 
Fig. 21 where the time evolution of the maximum value of 



FIG. 3: Hamiltonian constraint at t = 10 for simulations of 
a pure gauge wave with three different resolutions Ar = 0.1, 
Ar = 0.05, and Ar = 0.025 rescaled by the factors corre- 
sponding to second-order convergence. 
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FIG. 4: Time evolution of the maximum value of the radial 
shift /3 r (upper panel), and of the mass of the AH (lower 
panel) in the single puncture BH simulation. 



the radial shift j3 r is displayed in the upper panel. After 
an initial phase in which the maximum value of the shift 
vector grows rapidly, it settles to a value of ~ 0.15 and 
we find almost no drift until the end of the simulation 
at t = 2500. In the lower panel of Fig. [4j we show the 
time evolution of the mass of the apparent horizon (AH), 
defined as Mah = ^/A/Wtt, where A is the area of the 
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FIG. 5: L2-norm of the Hamiltonian constraint in the single 
puncture BH simulation. The insets show the L2-norm during 
the initial phase, and the L2-norm computed outside the AH 
respectively. 



AH. We notice that Mah is conserved well during the 
evolution and the error at t — 2500 is less than 0.2% 

The Hamiltonian constraint violation results are dis- 
played in Fig. [5j which shows the L2-norm of the Hamil- 
tonian constraint as a function of time. Both the ini- 
tial phase driven by the gauge dynamics (see upper inset 
in Fig. [SJ and the stationary phase are clearly visible. 
The lower inset in Fig. [5J shows that the L2-norm com- 
puted outside the AH is about three orders of magnitude 
smaller than the L2-norm computed in the whole grid, 
which is to be expected as the largest spatial violation of 
the constraint occurs near the puncture (due to the finite 
differencing of the irregular solution) . 



C. Spherical relativistic stars 

For our first numerical simulation of the coupling of 
Einstein equations and the general relativistic hydro- 
dynamic equations, we use the Tolman-Oppenheimer- 
Volkoff (TOV) solution. We focus on an initial TOV 
model that has been extensively investigated numerically 
by @, H3 ■ This model is a relativistic star with poly- 
tropic index N = 1, polytropic constant k = 100 and cen- 
tral rest-mass density p c = 1.28 x 10~ 3 , so that its gravi- 
tational mass is M = 1.4, its baryon rest-mass M» = 1.5 
and its radius R — 9.59. 

We evolve these initial data with our non-linear code 
until t — 3000 (~17 ms). In the upper panel of Fig. [6] we 
plot the time evolution of the central rest-mass density 
for a simulation with Ar = 0.025 and N r — 4000 until 
t = 1000. In the inset we show the same quantity for the 



FIG. 6: Upper panel shows the time evolution of the nor- 
malized central density for an M = 1.4, k = 100, N = 1 
polytrope. Power spectrum of the evolution of the central 
rest-mass density is shown in the lower panel. F, HI and H2 
represent the frequency of the fundamental mode and the first 
two overtones computed by [241 ]. 



whole evolution. We observe that the truncation errors 
at this resolution are enough to excite small periodic ra- 
dial oscillations, visible in this plot as periodic variations 
of the central density. We see that the damping of the 
periodic oscillations of the central rest-mass density is 
very small during the whole evolution, which highlights 
the low numerical viscosity of the implemented scheme. 

By computing the Fourier transform of the time evolu- 
tion of the central rest-mass density we obtain the power 
spectrum, which is shown with a solid line in the lower 
panel of Fig. [51 while the dashed vertical lines indicate the 
fundamental frequency and the first two overtones com- 
puted by [13]. Note that the locations of the frequency 
peaks for the fundamental mode and the two overtones 
are in very good agreement, the relative error in the fun- 
damental frequencies being less than 0.1%. 

The result of this simulation shows the ability of the 
scheme to maintain the numerical stability in long-term 
non-vacuum regular spacetime simulations in spherical 
coordinates without the need of an additional regulariza- 
tion at the origin. 



D. Gravitational collapse of a marginally stable 
spherical relativistic star 

We next test the capability of the code to follow BH 
formation with the gravitational collapse to a BH of a 
marginally stable spherical relativistic star. For this test, 
we consider a k = 100, JV = 1 polytropic star with 
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FIG. 7: Time evolution of the normalized central density (up- 
per panel) and mass of the AH in units of the ADM mass of 
the system (lower panel) for the collapse of a marginally stable 
spherical star to a BH. 
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FIG. 8: L2-norm of the Hamiltonian constraint in the collapse 
simulation. The vertical dashed line indicates the time of AH 
formation. The inset shows the L2-norm computed outside 
the AH. 



central rest-mass density p c — 3.15 x 10~ 3 , so that its 
gravitational mass is M = 1.64 and its baryon rest-mass 
M* = 1.79. In order to induce the collapse of the star, 
we initially increase the rest-mass density by 0.5%. 

We present numerical results for a simulation of the 
gravitational collapse of a marginally stable spherical rel- 



ativistic star performed with resolution of Ar = 0.125. 
We use the gauge conditions given in Eqs. (|2.20l) - (|2.22l) . 
We plot in Fig. [7] the time evolution of the normalized 
central density until t = 300 (upper panel), and of the 
mass of the AH in units of the ADM mass of the sys- 
tem until t = 500 when we stopped the simulation (lower 
panel). Overall, as the collapse proceeds the star in- 
creases its compactness, reflected in the increase of the 
central density as shown in the upper panel. The most 
unambiguous signature of the formation of a BH during 
the simulation is the formation of an AH. Once an AH 
is found by the AH finder, we monitor the evolution of 
the AH area, and also of its mass which is plotted, in 
the lower panel of Fig. [7] This panel shows that approx- 
imately at t ~ 167, an AH is first found and that the 
mass of the AH relaxes to the ADM mass of the system. 
The difference in the ADM mass and the mass of the AH 
at t = 500, is about 0.2%. 

In Fig[5]we plot the L2-norm of the Hamiltonian con- 
straint in the collapse simulation. The vertical dashed 
line indicates the time of AH formation, and the in- 
set shows the L2-norm computed outside the AH. The 
largest violation of the constraint occurs during the AH 
formation, and afterwards the value of the L2-norm set- 
tles to ~ 10~ 3 . As in the case of a Schwarzschild BH, the 
L2-norm of the Hamiltonian constraint computed outside 
the AH is about two orders of magnitude smaller than 
the L2-norm computed in the whole grid. 

Results of this simulation indicate that the numeri- 
cal scheme to integrate the evolution equations in time 
can handle accurately the transition between a regular 
spacetime (that of the star) and a irregular spacetime 
containing a puncture singularity at r = 0. 



VI. CONCLUSIONS 

In this paper we have presented a numerical code solv- 
ing the BSSN equations in spherical symmetry and the 
general relativistic hydrodynamic equations written in 
flux-conservative form. A key feature of the code is that it 
uses a second-order PIRK method to integrate the evolu- 
tion equations in time. This numerical scheme has proved 
to be crucial and sufficient to obtain the desired stabil- 
ity without the need for a regularization scheme at the 
origin. 

We have performed and discussed a number of tests to 
assess the accuracy and expected convergence of the code, 
namely a pure gauge wave, the evolution of a single BH, 
the evolution of spherical relativistic stars in equilibrium, 
and the gravitational collapse of a spherical relativistic 
star leading to the formation of a BH. We remark that, 
to our knowledge, we have presented the first success- 
ful numerical simulations of regular spacetimes (vacuum 
and non-vacuum) using the covariant BSSN formalism 
in spherical coordinates without the need for a regular- 
ization algorithm at the origin (or without performing a 
spherical reduction of the equations [28l. |29|]). 
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In addition, curvilinear coordinate systems facilitate 
the use of non- uniform radial grids (i.e., logarithmic ra- 
dial coordinate) to achieve the required high resolution 
near the origin while still keeping the outer boundaries 
sufficiently far away. This is particularly useful if one 
aims to study astrophysical phenomena like the gravita- 
tional collapse or the dynamics of accretion disks around 
BHs. Such approach is simpler, and likely computation- 
ally less expensive, than the adaptive mesh refinement 
techniques used in 3D codes in Cartesian coordinates. 

We note that, unlike with the Fully Constrained For- 
mulation [bl ] in which some of the equations take an 
elliptic form and where a similar PIRK has been suc- 
cessfully tested, the BSSN formulation is purely hyper- 
bolic, and yet the application of the PIRK method has 
proved very robust and provided the numerical stabil- 
ity necessary to perform long-term simulations of regular 
spacetimes in curvilinear coordinates. The work we have 
presented also paves the way for future comparisons of 
the performance in curvilinear coordinates between the 
BSSN and the Fully Constrained Formulation system. 
Moreover, the application of the PIRK method to the 
BSSN equations in 3D in such coordinate systems should 
be rather straight forward, and we aim to investigate this 
in a future work. 



the source terms included in the explicit or partially im- 
plicit operators are detailed. 

Firstly, the hydrodynamic conserved quantities, a, b, 
X, a and /3 r , are evolved explicitly, i.e., all the source 
terms of the evolution equations of these variables are 
included in the L\ operator of the second-order PIRK 
method. 

Secondly, A a and K are evolved partially implicitly, 
using updated values of a, a and b; more specifically, 
the corresponding L 2 and L3 operators associated to the 
evolution equations for A a and K are: 



L 



2(A a 



= - ( V r V r a - iv 2 a 



n ( K-^R 



L 3(A a ) = P r d r A a + aKA a - 16ira(S a - S b ), 
L 2 (K) = -V 2 a, 

L m) = p r d r K + a(A 2 a + 2Al + ^K 2 ) 
+ 4ira(E + S a + 2S b ). 



(Al) 

(A2) 
(A3) 

(A4) 



Then, A r is evolved partially implicitly, using updated 
values of a, a, 6, f3 r , conformal factor, A a and K; more 
specifically, the corresponding L 2 and L3 operators asso- 
ciated to the evolution equation for A r are: 
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Appendix A: Detailed source terms included in the 
PIRK operators for the evolution equations 



L 



2(A r 



L 



1 



d 2 r /3 r 



M- 

b \ r 



3o, 



2 Aa 
-(A a d r a + ad r A a ) - — (A a ~ A b ) 
a rb 



a 



d r A a - -^d r K + 6A a d r x 



+{A a - A b - + — 
1 r 



2(7 



3(Ar) = p r d r A r - A r d r (3 r + yA r V m ^ 
+ 2aA a A r -8Trj r — . 



(A5) 



(A6) 



The evolution Eqs. ([221), (|2T9]> , (I2~T1) . (|2~1B . 

(|2~T7|) . (PZ2H)) - ([2T2"2"]) , are evolved using a second-order 
PIRK method, described in Sec. III. In this Appendix 



Finally, B r is evolved partially implicitly, using up- 

3 - 

dated values of A r , i.e., L 2 (B r ) — ~7^tA r and L^ Br ) = 0. 
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